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We present results of simulations of a quantum Boltzmann 
master equation (QBME) describing the kinetics of a dilute 
Bose gas confmed in a trapping potential in the regime of Bose 
condensation. The QBME is the simplest version of a quan- 
tum kinetic master equations derived in previous work. We 
consider two cases of trapping potentials: a 3D square well 
potential with periòdic boundary conditions, and an isotropic 
harmònic oscillator. We discuss the stationary solutions and 
relaxation to equilibrium. In particular, we calculate particle 
distribution functions, fluctuations in the occupation num- 
bers, the time between collisions, and the mean occupation 
numbers of the one-particle states in the regime of onset of 
Bose condensation. 



I. INTRODUCTION 



In the previous two papers [[lj|2|], which we will refer 
to as QKI and QKII, a fully quantum mechanical ki- 
netic theory for Bose gases was developed. One of the 
simplest versions of the quantum kinetic master equa- 
tion (QKME) neglects all spatial dependence, and yields 
a master equation, which we have named the quan- 
tum Boltzmann master equation (QBME). In contrast to 
the familiar quantum Boltzmann equation (QBE) 
which is an equation of the single particle distribution 
function, the QBME is an iV-atom stochastic equation. 
The aim of the present paper is to present results of nu- 
merical simulations of this equation for finite size systems 
consisting typically of a few hundred àtoms. Although 
the exclusion of the spatial dependence is an extreme 
simplification, these simulations will give us a first orien- 
tation about the kind of solutions the QKME will yield. 
These simulations can thus serve as a guideline for the 
type of approximations of the QKME one may use to 
find numerical solutions of this much more interesting, 
but unwieldy, equation. 

Furthermore, we will concentrate our attention on 
those results of the QBME which cannot be obtaincd 
using equations like the quantum Boltzmann equation 
(QBE) . We also restrict our work to the region of tem- 
peratures which are less than or not much higher than 
the critical temperature of the gas, because at much 
higher temperatures quantum effects do not play a crucial 
role and simulations of the classical Boltzmann equation, 
which is vàlid in that case, have already been performed 
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The QBME is a genuine iV-atom equation like the 
QKME, but it neglects all the coherences contained in 
the QKME — it is thus intermediate between the QKME 



and the description of the system with kinetic equations 
for single-particle distribution functions. Its irreversibil- 
ity comes from the Markov assumption employed in de- 
riving the QKME. 

The paper is organized as follows. In Sec. ÏJ we re- 
vicw the derivation of the QBME in QKI JÏJ, discuss 
properties of the QBME, and compare it with the QBE. 
Furthermore, we give a brief description of the simula- 
tion algorithm. In Secs. [II and IV we apply the QBME 



to study a Bose gas confmed in a 3D box and in a 3D 
harmònic oscillator. In particular, we discuss simulations 
results for thermodynamic quantities, the mean time be- 
tween collisions, and the fluctuations of the occupation 
numbers of the condensate. For the 3D harmònic oscil- 
lator we also simulate a gas that is evaporatively cooled. 



II. THE QUANTUM BOLTZMANN MASTER 
EQUATION 

In this section we will first summarize the derivation of 
the quantum Boltzmann master equation as given in QKI 
0. Furthermore, we discuss properties of this equation 
and its solutions which are relevant for our numerical 
studies presented in Secs. [II and and conclude with 
a comparison of the QBME with the QBE. 



A. Derivation and validity of the QBME 

The second quantized form of the Hamilton operator 
for a Bose gas with pair particle interaction can be writ- 
ten H = Hq + Hj , where 



Ha 



2J fuJx, 



Hj 



E 

mi,m2,m 3 ,m4 



in i .m; .rri3 ,m4 



a mi a m 2 a m 3 O r 



(D 



(2) 



Here Hq is the system Hamilton operator of the non in- 
teracting Bose gas where is the creation operator of 
a particle in the eigenstate of Hq labeled ni; with energy 
Tujj mi - The trapping potential is included in Hq. 

The interaction Hamiltonian Hi describes two body 
interactions in the Bose gas. In the regime we want to 
study, only s-wave scattering plays an important role, 
allowing us to write 
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In Eq. (||) ^f mi (x) denotes an eigenfunction of the trap- 
ping potential in coordinate space labeled by quantum 
numbers i. Below we will specify the potential to be a 
3D box with periòdic boundary conditions (Sec. [II ), or a 
3D isotropic oscillator (Sec. IV), and will give expressions 
for the matrix elements U, 



mi ,1112 ,1113 ,1114 



for these specific 
cases. For convenience we will use the notation i instead 
of rrii below. The scattering length of the gas is a and the 
mass of the gas partides is m. We will treat systems with 
finite number of partides N. This is the starting point 
from which the QKME is derived in QKI. The following 
assumptions and approximations are made: 



1. The forward scattering temis 

All the terms of Hj (see QKI Eq. (67)) which commuto 
with the system Hamiltonian H describe forward scat- 
tering and give rise to the mean field. These terms can 
be included with Hq. Forward scattering does not change 
the occupation of the one-particle eigenstates, so we will 
neglect the influence of these terms on the eigenstates of 
H in the simulations. 



2. The collision terms 

The remaining terms in Hj describe col·lisions which 
change the occupation numbers of the one-particle states 
of the trapping potential. We assume that this part of 
Hj can be treated perturbatively, using the Born approx- 
imation and the Markov approximation (QKI Sec. IV C 
3). The Born approximation is vàlid when the interac- 
tion between the partides is small compared to the sys- 
tem Hamiltonian Hq || . For the Markov approximation 
to be vàlid it is required that the frequency spectrum is 
effectively continuous which means that the separation 
between the energy levels is much smaller than the en- 
ergy range of occupied states. The use of the Markov 
approximation gives the QKME its irreversible charac- 
ter. We will neglect the influence of collisional shifts on 
eigenstates of H. 



3. Reduction of the QKME to the QBME 

To reduce the QKME to the QBME it is assumed 
that the coherent terms (i.e. Hamiltonian terms in QKI 
Eq. (77)) can be neglected. The QBME is an equation 
for the diagonal elements w n = (n|p|n) of the density 
operator and takes the following form (QKI Eq. (101)): 
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X {mn.2 (n 3 + 1) (n 4 + 1) [lV n - W n+ei23i ] 

+ (ni + 1) (n 2 + 1) n 3 n 4 [w n - w n -e 1234 ]} • (4) 

Here |n) = \iïq, ni, ri2, ...) is a Fock state of the iV-particle 
system, giving the occupation numbers rn of the eigen- 
states í'i(x) and n denotes the vector consisting of the 
occupation numbers rij. The vector ei2 3 4 is defined sim- 
ilarly to n as 

ei234 = [0,...0,Í,0,...0,Ï,0...0,-1,0,...0,-1,0,...0], (5) 

which describes two particle col·lisions. The state |n — 61234) 
can thus be reachcd from |n) by the collision 1+2 — > 3+4. 

The (5-function in the discrete sum of the QBME ([|) 
has its origin in the use of the Markov approximation as 
outlincd in QKI. Since we do not replace these sums by 
integrals in our simulations this <5-function requires inter- 
pretation. We concentrate energy regions of Ae to one 
single discrete energy level. The energy interval is de- 
scribed by the properties of the closest one-particle state 
of Hq. The choice of Ae depends on the trapping poten- 
tial and is such that each of the one-particle states serves 
as one of the discrete energy levels with energy ei that 
determine the properties of partides within the energy 
range \e% — Ae/2, + Ae/2]. Implicitly this includes the 
interpretation of rij as being an integral over a smooth 
distribution function /(e): 



i- Ae/2 Ae 



(6) 



where /(e) gives the number of partides occupying a 
state with energy e. Among the degenerate one-particle 
eigenstates the partides are distributed according to sim- 
ilar arguments as in Eq. (||). The (5-function in the 
QBME (|j) has, therefore, to be interpreted as 



Ae 



(7) 



and w n in the QBME (||) is the probability of finding m 
partides within the energy interval [e.; — Ae/2, ej + Ae/2]. 
5 XjV denotes the Kronecker delta. 

In QKI it is shown that the kernel of the integral where 
the Markov approximation is made (QKI Eq. (68)) has 
a width given by the temperature h/kT. This width 
determines the range of possible outeomes of a 

collision. As long as h/kT is much smaller than the 
time between two col·lisions the free evolution after the 
kernel has reached zero will fix the energy of the partides 
within a range of %/t co \\ before the next collision oceurs. 
We already assumed that it is possible to describe the 
system in terms of one-particle eigenstates of the trap- 
ping potential which is only vàlid if the level broadening 
coming from the collisions is much less than the level 
spacing. Hence, we are able to decide which of our one- 
particle states describes the properties of a particle best 
before this particle collides again. 
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4- Discussion 

Since there are no classical assumptions in deriving 
the QBME (Q), it should be vàlid cvcn when thc Bosc 
gas becomes degenerate, within the limits of the approx- 
imations made in its derivation. The quantum statistics 
is contained in the 1 + rii factors in Eq. (|j) . This allows 
us to study the onset of BEC, in the sense of obtaining 
a macroscopic occupation in the ground state [Q, and 
in particular finite number effects which are important 
when the number of àtoms is not large. 

The QBME (|) is a full A-particle equation in the form 
of a stochastic master equation which describes N parti- 
des interacting with each other by two particle col·lisions. 
These col·lisions are responsible for the equilibration pro- 
cess. In contrast, the QBE (see Sec. |II Dj ) considers the 
motion of one particle interacting with a mcan distribu- 
tion of the other àtoms in the gas. 

No mcan ficld cffccts arc included in the present form 
of the QBME (^). As soon as the temperature T of the 
gas is far below the critical temperature T c and most of 
the partides have aceumulated in the ground state, the 
mean field produced by these condensed partides must 
be taken into account. Furthermore, the derivation of 
the QBME assumes that the width of the energy levels 
and collisional shift in addition to the mean field is small 
relative to the level spacing Ae. 



B. Quantities of interest 

For comparison with the simulations discussed in the 
following sections, we summarize below properties of the 
stationary solutions, the particle distributions and colli- 
sion times. 



1. Stationary solution 

The QBME conserves energy E and number of parti- 
des N. According to QKI the stationary solution of the 
QBME (|) is 



constant, 



(8) 



corresponding to a microcanonical ensemble. 

We will also compare our simulations results with the 
grand canonical ensemble. For the mean occupation 
numbers one obtains (compare QKI Sec. V A 2) 



(m) = 



1 



exp 
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(9) 



In this case T is the temperature and /i is the chemical 
potential of the system in the grand canonical ensemble. 
Given the mean energy of the system E and the mean 
number of partides N we can solve the two equations 



N = J2 
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for fi and T numerically. We will compare this result 
below with the one we get from our simulations. In thc 
framework of the QBME these grand canonical results 
are obtained if we assume that in steady state the expec- 
tation vàlues of the rij factorize (which is an approxima- 
tion). 



2. Particle distributions 

The QBME is a stochastic equation for the diagonal 
elements of the density operator in the basis of the eigen- 
states of H . We are interested in calculating the proba- 
bility distribution of partides in the one-particle states. 
They are defined as 



(ii) 



and give the probability of finding j partides in the one- 
particle eigenstate labeled i. The sum runs over all n 
with = N and h Y]j uijUi = E, the constant num- 

ber of partides in the gas and the energy of the system, 
respectively. We will compute these distributions for the 
3D box in Sec. [III B 4 . 

For highly excited states i, whose mean occupation 
number is much less than 1 the probability Wi(j) is only 
substantially different from zero for j = and j = 1, 



which leads to (nf ) 



On the other hand, 



approximate expressions can be derived for low lying 
states, including the ground state, on the following argu- 
ments. Assuming that there is no restriction on how the 
partides are distributed among degenerate energy levels 
we can write Wi(j) in terms of energy levels M 



Wj(j) 
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(12) 



Here ï are sets of indices of degenerate eigenlevels, gj 
is the number of elements of l, and grni gives the total 
number of partides in the states l € l. The normalization 
constant is denoted by Z ', and n is a vector containing the 
nj. This formula is only approximate because it includes 
configurations of the system that cannot oceur in the sim- 
ulations since they are not connected by collisions with 
the initial configuration. For small temperatures the sum 
in Eq. (|ï^) is readily calculated numerically, an d we wi ll 
compare this with our simulation results in Sec. III B 4. 
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3. Collision time 

For any given configuration n of the system we calcu- 
late the sum over all the transition matrix elements for 
collisions that can occur. This sum is the value of the 
right hand side of Eq. (|]) for a given n, the correspond- 
ing w a = 1 and all the other w m equal to zero. A single 
possible collision 1+2^3+4 contributes to this sum 



P(12 — ^ 34) 



47T 

hAe 



Iü'umI nina (n 3 + 1) (n 4 + 1) . (13) 



where the factor of 4 is due to different permutations of 
the indices which describe the same collision. We call 
a collision possible if it conserves energy and C/1234 7^ 0, 
and we call P(12 — > 34) the transition probability per unit 
time for this particular collision. 



C. Ergodic approximation 

In solving the QBE or the classical Boltzmann equation 
it is a common approximation to simplify this equation 
by an ergodic assumption ||MEBpj . In a classical con- 
text this corresponds to the assumption that the phase 
space density / t (p,x) only depends on the energy e of 
the partides at position x with momcntum p at time t. 
Quantum mechanically, it is postulated that degenerate 
energy levels carry equal populations at all times, i.e. the 
populations of degenerate eigenlevels equalizc on a time 
scale much faster than collisions between levels of differ- 
ent energies. This implies that the occupation numbers 
m in the QBME should be replaced by 



ffï^ 



(14) 



Hcrc wc dcfinc sets of indices i that contain all the in- 
dices of one partides states with the same energy hiúi\ 
and g- t is the degeneracy factor of states with energy lu^. 
We note that the n\ are no longer integers. In our sim- 
ulation this corresponds to a distribution function which 
is completely specified by the occupation numbers of (the 
block of) degenerate energy levels, i.e. w n — > u> fi where 
n is a vector containing the number of partides in the 
degenerate eigenlevels nj. Removing or adding a parti- 
cle to a state i changes nj by l/gj. Therefore, we use a 
vector ej234 which is defined by 

e I53 4 = [0,...0, 1/51,0,.. .0,1/52, 0... 

...0,-1/53,0, ...0,-1/54, 0,-0], (15) 

to describe collisions in the ergodic case. 

Using these definitions we can write the ergodic form 
of Eq. (|4|) in the following way 



1234 



/ 



s {h,{u)i + Lü 2 - lü 3 - UJ4)) 1 u 123 4 r 



iei,2e2 

3E3.4S4 



x {nin2 (ng + 1) (ni + 1) [w & - Wn+ e - l3 - 3l ] 

+ { n ï + 1) («2 + 1) ™3™4 [w & ~ %-ej K ] } ■ (16) 

Transition probabilities P(Ï2 — > 34) are calculated ac- 
cording to 



P(12 34) = P ( 12 34 ) 



(17) 



iei,2e2 

3£3,4e4 



where the sum runs over all the elements of a particular 
set of degenerate states. Collisions which do not change 
the energy distribution are thus no longer taken into ac- 
count. 

We note that the ergodic assumption yields the correct 
steady state distribution, but we expect differences in the 
details of the dynamics. A comparison of the kinctics 
with and without the ergodic assumption will be given 
in the case of a 3D box in Sec. III; our simulation results 



IV will be based on 



for the harmònic oscillator in Sec. 
the ergodic approximation. 



D. Comparison between the QBME and the QBE 

Recent work of kinetics in relation to Bose conden- 
sation in trapping potentials by Holland and collabora- 
tors (5| is based on the QBE with an ergodic assumption 
(for a classical Boltzmann equation see also ||). The 
derivation of the QBE is based on factorizing mean vàl- 
ues (nin 2 ...ni) = (m) (ra 2 )...(raj) with (m) = J2 n n i w ^- 
In the ergodic approximation one obtains 

47T 

9ï(n%) = — g{-(ni)(n5)«ng) + l)«n a ) + l) + 

234 

((n I ) + l)((n 5 ) + l)(ng)(n 4 )} 



|t r i234| 2 <5 (%(W1 +ÜJ 2 - Ld 3 - LO4)) | . (18) 

leI,2G2 
V3£3,4e4 

In this equation (nj) is the mean occupation of the de- 
generate states. The discretc (nj) replace the particle 
distribution function /(e) used in the classical version 
H . The QBE describes the time evolution of the single 
particle distribution function in the mean distribution 
of the other partides. In contrast to simulations of the 
QBME, fluctuations in the occupation numbers are thus 
not described by the QBE. Moreover, it is not possible 
to simulate systems far from equilibrium where the fac- 
torization of the mean vàlues is not vàlid. On the other 
hand, the QBE has the advantage that it allows simula- 
tions with much larger particle numbers than the QBME. 
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E. Simulation of the QBME 

Since the QBME is a stochastic master equation for 
(quantum mcchanical) occupation probabilities, we can 
simulate its time evolution as a series of jumps. A jump 
describes the collision of two partides 1+2 — » 3+4, which 
is represented as an instantaneous change of the corre- 
sponding occupation numbers. The simulation method 
is often used for rate equations and works as follows: 

1. Take an initial configuration of partides, n, (rep- 
resenting an initial density operator pg(t = 0) = 
|n)(n|), where the energy E and the total number 
of partides N are fixed. 

2. Calculate all the transition probabilities per unit 
time P(12 — > 34) for the given n. 

3. The total collision rate is now proportional to the 
sum over all transition probabilities per unit time. 

4. The next jump occurs at time t m since the last 
jump, which can be calculated by choosing a ran- 
dom number r €]0, 1] from a uniform distribution 
and using 



ln(r) 



Ei234 ^(12^34) ■ 



(19) 



5. All the possible collisions are lincd up with the 
length P(12 — > 34). Another random number 
s s]0, X)i234 ·f(12 34)] is chosen from a uniform 
distribution. The transition selected by this ran- 
dom number s gives the particular collision 12 — > 34 
which occurs. 

6. The last step now is to set t := t + t m , |n) := 
|n - ei234) and p e := |n)(n|. 

7. Go back to 2. 

8. Repeat this simulation to obtain p = constant x 
J2e Pi- 
hi every collision only four of the occupation numbers 

are changed and therefore only fcw of the transition ma- 
trix elements are modified by the change in the occupa- 
tion numbers. Thus it is not necessary to calculate all 
transition probabilities after each step, since only those 
involving the m, ri2, «3, n<i which define the collision 
will have been changed. (This is, however, more compli- 
cated than is the case for the Boltzmann master equation, 
where the 1 + íij factors do not occur) 

For integer occupation numbers it is not possible to 
neglect the (1 + n,) factors above a certain energy by 
arguing that the mean occupation of highly excited states 
is much smaller than one. For highly excited states these 
factors are either 1 or 2, etc. and cannot be replaced by 1. 
We have to account for them regardless of the energy of 



the one-particle states involved into the collision. In our 
simulation method we do not restrict the number states 
available for the partides of the gas. We keep track of 
each of the partides rather than of a certain number of 
one-particle states. This limits the number of partides 
we are able to consider. 



III. 3D SQUARE WELL POTENTIAL WITH 
PERIÒDIC BOUNDARY CONDITIONS 

A. Description of the system 

First, we will simulate the QBME for a 3D cube of 
length L with periòdic boundary conditions. This cor- 
responds to the simplest version of the QBME. In the 
language of QKI Sx — L is the length of the phase space 
celis. From this we immediately find the spacing of the 
celis in momentum to be Sp = 2nh/L, which is equal 
to the momentum spacing of the discrete energy levels 
in the box. The wavelet functions (introduced in QKI 
Eq. (26)) are therefore reduced to 



zkx 



Vk(x) = 



L 3 ' 2 



(20) 



We have dropped r of QKI in the equation above because 
there is only one phase space celi in coordinate space. 
The wave numbers k take on the discrete vàlues 



k = — m 

L 



(21) 



where m is a vector consisting of integer vàlues. Since 
our system has the finite volume L 3 the wavelet functions 
are orthogonal in the following sense 



J d 3 x w ki (x) v£ ( (x) = Si.i 



(22) 



With these wave functions we can now calculate t/1234 to 
be 



Pl234 = 



ATrh 2 a 
mL 3 



J mi+m2,m 3 +m4 • 



(23) 



In the case of the 3D box Ae = (2h 2 1: 2 ) / (mL 2 ) . Us- 
ing a = 87ra 2 pi for the cross section, n — N/L 3 and 
Vi = (2irh)/(mL) which is the magnitude of velocity of 
a particle in the first excited state we get 



P(12 -> 34) = anvi 



-'mi +rri2 ,«13 +m4 



Nit 

nin 2 (n 3 + 1) (n 4 + 1) . 



(24) 



The number of possible collisions is restricted by two 
Kronecker delta functions that ensure energy and mo- 
mentum conservation. The overlap integral is P1234 = 
{Airh 2 a) / (mL 3 ) for all the possible collisions, and does 
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not depend on energy or momentum of the involved one- 
particle states. 

In semi-classical treatments of the QBE, the ergodic 
assumption is often made |||4j . The density of states is 
approximated to be proportional to y/e. It is then shown 
that the transition matrix elements are proportional to 
y/e m i n , where e m i n is the minimum energy of the colliding 
partides (compare Appendix |a|). In case of a smoothly 
varying, strictly decreasing hmction /(e), one can there- 
fore argue that most of the col·lisions happen between 
partides with almost the same energy. In the cases we 
are interested in we cannot make these assumptions. The 
occupation numbers can vary strongly and the degener- 
acy of states which we count exactly is not proportional 
to ye in the energy range in which our simulations are 
performed. 

B. Results of simulations 

All the simulations we report contain a statistical error. 
Unless this statistical error is given explicitly it is less 
than 5%. 



The degenerades gi in Eq. ( |ï(ï| ) are calculated by count- 
ing all ways of combining different integer numbers mf , 
to ,m| consistent with the defmite energy 

fi** = («) 2 + K) 2 + K) 2 ) • (25) 

The simulation and the grand canonical result both 
give a higher number of partides in the condensate than 
expected from the thermodynamic result. The results for 
finite number of partides, however, approach the thermo- 
dynamic limit with increasing N very quickly. Around 
the critical temperature there is a slight deviation of the 
simulated results from the grand canonical results | {Ï3| , 
whereas for T <^T C there is almost no difference. This is 
due to a bigger statistical error in the simulation because 
of the large fluctuations in the region around the criti- 
cal temperature. Note also that we are comparing two 
different statistical ensembles, and that for finite systems 
we would not expect exact agreement between the results 
from different ensembles. 



2. Occupation of the ground state 



1. Thermodynamic quantities 

There are two ways of computing the stationary solu- 
tion of the QBME (||) . The first is to calculate it directly 
from Eq. (j8|); this is only feasible for very few àtoms. 
The second possibihty is to obtain the stationary solution 
from simulations by assuming that the time average over 
a sufficiently long time period equals the ensemble aver- 
age. To find this time we wait until the simulation results 
agree with a Bose-Einstein distribution Eq. (||) . This also 
allows us to assign a temperature T to the systcm. All 
the results are sealed to the critical temperature in the 
thermodynamic limit T c = 2h 2 n/(mL 2 k)(N/((3/2)) 2 / 3 
|jïïf . There are three parameters of the system; E, L and 
N which give a certain T, T c and N in thermodynamic 
equilibrium. In the simulation we fix E and N, and the 
sealing to the critical temperature is equivalent to sealing 
to a certain particle density N / L 3 in the box. 

The expression T c for the critical temperature is, of 
course, only vàlid in the thermodynamic limit because 
in the derivation |ï^ ] sums over energies are replaced by 
integrals which over- or underestimate the sums for finite 
systems depending on the density of states. In the ther- 
modynamic limit the energy spectrum becomes continu- 
ous and summing yields the same result as integrating. 
For a finite number of partides we therefore do not ex- 
pect the critical temperature and the condensate fraction 
vs. temperature to be the same as in thermodynamic 
limit. 

Fig. [ï] shows the comparison of the results from the 
simulations and the grand canonical expression Eq. (10). 



We want to investigate the sealing of the one-particle 
state occupation with the number of partides in the gas. 
For BEC we expect || the occupation of the ground state 
to be 




3/2N 



while for excited states 



(ni) T mk 
1 3 ~-~L X 2n 2 h 2 ia 2 



(26) 



(27) 



Fig. g shows the occupation numbers of the ground state 
and the first excited state. At a given T/T c the number 
of partides in the ground state increases linearly with 
the total number, whereas the slope of the occupation 
of the first excited state becomes smaller with increasing 
number of partides. From this numerical result we con- 
clude that the QBME does really describe a macroscopic 
occupation and is consistent with expecting BEC below 
the critical temperature T c . 



3. Collision times with and without the ergodic assumption 



O bvious ly, the results given in previous Secs. III Bl 
and III B 2 are the same with or without use of the ergodic 
assumption. This is expected, because in thermal equi- 
librium all degenerate eigenlevels should have the same 
occupation number even without the ergodic assumption. 
The mean time between two collisions in thermal equi- 
librium is computed by taking the average over all the 
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calculated times t m from Eq. (^9|). This time has to be 
multiplied by N/2 because one particular particle is in- 
volved in one out of N/2 col·lisions. 

In Fig. D we plot the mean collision time í"® n for onc 
particle versus the temperature without the ergodic as- 
sumption . The classical elàstic mean collision time cal- 
culated from í£ oll = (crol'T) , where vt is the mean ther- 
mal velocity of the gas vt — N^ 1 J^. -J (2hu>i) / (m) (rij) , 
and (rii) is the mean occupation of the i-th energy level 
obtained from the simulation. As soon as the gas be- 
comes degenerate the 1 + rij factors in Eq. (||) become 
important and increase the collision rate compared to the 
classical case. For temperatures close to zero the collision 
time increases again because there are only few partides 
outside the condensate which can take part in collisions. 
Fig. |] shows the comparison between the curves for the 
ergodic collision time í° oll and f"° u . For very small tem- 
peratures, the ergodic assumption allows for collisions 
which can not occur in the non-ergodic case because the 
corresponding states are not occupied. As soon as the 
temperature is close to T c , those collisions in the non- 
ergodic case that only change the direction of the mo- 
mentum of an individual particle decrease í^ n compared 
to the ergodic case. Since this type of collisions leaves the 
energy of the partides unchanged they are not included 
in the ergodic calculations. 



4- Particle distributions 

While the mean vàlues for the occupation numbers are 
easy to calculate it takes more effort to find the parti- 
cle distribution Wi(J) of the one-particle states. There 
are two ways of calculating these distributions. Either we 
calculate the time a state was occupied by a certain num- 
ber of partides (time average) or we record the numbcr 
of partides in that state after a certain time for many 
different trajectòries (ensemble average). Both results 
nccd not necessarily be the same unless the system has 
the stationary solution Eq. (||). We used both methods 
to calculate particle distributions for the condensate and 
some of the excited states for different temperatures. 

In Fig. U particle distributions for the ground state are 
plotted. Particle distributions of the condensate are well 
approximated by a Gaussian for temperatures below T c . 
However, they are not complctdy symmetric around the 
mean value like the Gaussian there is a slight asymmetry 
which increases with temperature. The shape of the dis- 
tribution changes close to the critical temperature. For 
N = 500 at T = 1.1T C the distribution develops a second 
local maximum at N c = and at T = 1.2T C the peak 
at finite number of condensate partides has disappeared. 
Well above T c at T — 1.7T C it agrees with a Bose-Einstein 
distribution 



The particle distribution of the first excited state in 
Fig. D can be approximated by the Bose-Einstein distri- 
bution (53) for r«T c and T > T c . Particle distribu- 
tions of highly excited states agree with the Bose-Einstein 
probability distribution ( p8| ) at all temperatures. 

In Fig. Q we plot the Standard deviation <t(N c ) of 
the particle distribution of the condensate in thermal 
equilibrium. The error bars are calculated according 
to ^/(t((N c — (N c )) 2 ) / (N c ) which is the variance of the 
Standard deviation normalized to the mean number of 
partides in the condensate. This gives the mean devia- 
tion of the Standard deviation from its calculated value. 
For small temperatures a(N c ) rises almost linearly with 
temperature. The number of possible states with differ- 
ent number of partides in the condensate increases which 
leads to a larger width of the distribution. Close to T c we 
get a very broad particle distribution with a very large 
Standard deviation. For T > T c the Standard deviation 
will tend to go to the mean number of partides in the 
condensate which agrees with the fit to the Bose-Einstein 
distribution which has a Standard deviation of N C (N C + Í) 
going to iV c for N c -C 1. 

This also agrees with the calculation performed in 
Sec. [I B 2 for the case of very small mean occupation 



p(N c ) = (1 - v )if' 
with (JV C ) = 77/(1 - 77). 



(28) 



of a state. At small temperatures we calculate the parti- 
cle distribution in the condensate according to Eq. (IÏ2]). 
To compute the sum in Eq. ( |Ï2| ) we assume that most of 
the fluctuations come from exchange of partides of the 
condensate with the first few excited states. The parti- 
des in higher excited states should not have a significant 
influence on the fluctuations in the condensate, but they 
should ensure that partides in the lowest lying states 
can be distributed among degenerate eigenlevels, with- 
out restrictions due to conservation laws. The derivation 
of Eq. (jï^) is based on the assumption that the parti- 
cle distribution among degenerate eigenlevels is not re- 
stricted by conservation laws. Using Eq. (|ï^) we obtain 
particle number fluctuations of the condensate due to ex- 
change with low lying levels. In particular we calculate 
Wç,(j) from Eq. ( |l2| ) by taking into account the first sev- 
enteen energy levels. The particle distributions we get 
agree well with the ones from the simulations. In Fig. |ï] 
the results of both calculation methods are compared for 
temperatures T < 0.5T C . The crosses correspond to the 
numerical calculations on Eq. (|ï^) and agree well with 
the simulation results. 



5. Growth of the condensate 

Here we want to investigate how the condensate builds 
up when the simulation is started in a non-equilibrium 
distribution. As the initial state we choose a Gaussian- 
like distribution: We first distribute the partides ran- 
domly into states with energies between that of the first 
excited state, and twice the mean energy, and then move 
partides to higher or lower energy states until the given 



7 



fixed energy E of the system is exactly reached. 

Whenever possible we avoid putting partides in the 
condensate at the beginning of the simulation. As can 
be seen in Fig. || the condensate growth is well fitted by 
N c (1 — exp(— t/r)) where N c is the number of partides 
in the condensate in thermal equilibrium and the time 
constant r is found by fitting this function to the simu- 
lation. This holds as long as the fraction of condensate 
partides in thermal equilibrium is not much less than 
one. 



the mean magnitude of velocity of a particle in the 
ground state of the oscillator, and ffi h is the mean parti- 
cl e density if all the partides are within a cube of length 
\J (hn) / (muj) . This is the semi-classical expression ob- 
tained in Appendix |A| According to Ref. [|Ï0| numerical 
calculation shows that the this expression is a good ap- 
proximation even for low lying energy levels. 

B. Results of simulations 



6. Time to reach an ergodic distribution 

While with the ergodic assumption all degenerate lev- 
els are equally occupied at all times, in the non-ergodic 
case col·lisions themselves are responsible for equalizing 
the occupations of degenerate levels. To check the relax- 
ation time for a distribution to become ergodic we dis- 
turb a system in thermal equilibrium by putting all par- 
tides with energy Ae into two of the first excited states 
(with opposite momentum so that the total momentum 
is unchanged). As is shown in Fig. || the particle dis- 
tribution comes to equilibrium in approximately 10í"° n . 
Col·lisions, therefore, transfer the occupation between de- 
generate levels at a time scale of the order of the mean 
collision time in the gas. We conclude that for the er- 
godic assumption to be vàlid strictly speaking it is only 
reasonable to look at quantities that are mean vàlues over 
sevcral collision times. 



IV. 3D ISOTROPIC HARMÒNIC OSCILLATOR 

A. Description of the system 

In this section we will study Bose partides trapped in 
an isotropic harmònic trap with trap frequency lj. The 
vector |n) now gives the occupations of the trap levels, 
and Í7i234 contains the spatial eigenfunctions of the har- 
mònic oscillator. For the low lying levels these integrals 
can be evaluated numerically but for highly excited states 
it is difficult to get reliable results for í/i234- Therefore, 
we will limit ourselves to u sing the ergodic form of the 



QBME as explained in Sec. II C. As is shown in pTj, the 



transition matrix elements of transitions which change 
the energy distribution function can be approximated by 

-- -- 4ir mauj^fr h 
P(12 — > 34) = j^- ^ 3 ff m in(l234) 

n ï n 2Í n 3 + l)( n 3 + 1) 

= n h av (2N)-'9^ n{ïm 
nïn2(n g + l)(nj + l) (29) 

Here gk = (j + + 2)/2 is the degeneracy factor of 
the j-th. eigenstate with energy jhuj, vq — \J (4huj) / (nm) 



1. Stationary Solutions 

To obtain the grand canonical stationary solutions for 
the 3D harmònic oscillator we have to replace gj by in 
the Eq. jïo|). The critical temperature for an ideal Bose 
gas in a 3D isotropic trap in the thermodynamic limit 
(i.e. when the sums over the discrete energy levels are re- 

1/3 

placed by integrals) is given by T c = (Tno)/k (N/((3)) 1 
|Ï2| . Our simulation results for the condensate fraction 
versus temperature are shown in Fig. [ï^. The continuum 
approximation increases the condensate fraction for finite 
number of partides compared to the simulation results. 
The reason for this is that the density of states rises much 
faster than for the 3D box. As in the case of the square 
well potential, the results for the microcanonical simu- 
lations and the grand canonical calculations agree very 
well. Comparing the two curves for N = 500 of Fig. |l| 
and Fig. 10 we find that the phase transition is more pro- 
nounced in the harmònic oscillator compared to the much 
smoother transition for the 3D box. This behavior can 
also be seen by plotting of the energy versus temperature 
in Fig. O. There is a visible change in the slope of the 
energy for the harmònic oscillator even for N = 500. For 
the ideal gas the heat capacity has a jump at the critical 
temperature in the thermodynamic limit in the harmònic 
oscillator whereas in case of the 3D box only the slope of 
the heat capacity is discontinuous at the critical temper- 
ature This makes clear that the thermodynamically 
cxpected differences in the condensation process between 
the harmònic oscillator and the free gas can also be seen 
in finite systems for small particle numbers. 



2. Collision times 



We will now compare the mean collision time obtained 
in our simulations íj?° n with the elàstic collision time de- 
(íï/ l (7i>Th) -1 · We determine wxh and 



fined as t^ ü 

with the assumption that the kinetic energy of the parti- 
des are equal potential energy equally. Then we find 
for the mean density and the thermal velocity üh = 

(3iV)/(47r) ((mLü 2 )/(E 3/2 )) 3/2 and v Th = y/E 1/2 /m rc- 

spectively, with E s = {l/NJ2i {huiY (»»)) • 

In Fig. ^| we plot the mean collision time versus tem- 
perature. For temperatures higher than the critical tem- 
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perature the simulation agrees well with the classical re- 
sult (dashed curve). For temperatures far below T c , the 
result of the simulation is approximately equal to the 
dotted curve which we obtained by the assumption that 
the size of the cloud is the ground state size and only the 
thcrmal vclocity varies with temperature. Around the 
critical temperature the size of the cloud shrinks faster 
than expected from the classical approximation. 

3. Evaporative cooling 

Currently BEC is achieved in experiments by evapora- 
tive cooling, i.e. by removing partides with a high energy 
from the trap (for a review see jïij ) . Elàstic collisions be- 
tween the partides thermalize the particle distribution 
which leads to a decreasing temperature. To simulate a 
Bose gas that is evaporatively cooled we cut off the trap 
at a certain energy level E b (t), with E b {t) a given func- 
tion of time. Each particle which is scattered into an 
energy level above E b {t) after a collision is considered as 
lost. In our simulations we start with Nq — 800 partides 
in the thermodynamic equilibrium at a temperature of 
T w 1.4T C . Then all partides with an energy larger than 
Eb(t = 0) = 65?iu; are removed. During the simulation 
we decrcase E b exponentially according to 

E b (t) = (E b (0) - £,)e~ 7í + E h (30) 

where Ei = 8Hu>. In Fig. [Ï3| the total number of par- 
tides in the gas N and the number of partides in the 
condensate N c are plotted as a function of time for dif- 
ferent parameters 7. First the partides in the highest 
energy levels are evaporated quickly. During the cooling 
process, the collision time decreases by an order of mag- 
nitude as shown in Fig. nevertheless the number of 
partides evaporated per unit time does not increase dur- 
ing the cooling process. The reason is that most of the 
collisions occur between partides with almost the same 
energy and thus many collisions are necessary to redis- 
tribute the partides when some of them are evaporated. 
If the collision rate did not increase so rapidly partides 
might be lost from the trap faster than evaporative cool- 
ing is possible. As soon as the condensate builds up the 
mean collision time increases again. This expected be- 
havior agrees qualitatively with Fig. |ï^. 

In order for the evaporative cooling to be efhcient it 
is important to quickly put as many partides as possible 
into the condensate. We therefore have calculated the 
size of the condensate divided by the time needed to reach 
90% of the equilibrium condensate fraction for different 
vàlues of 7 as shown in Fig. |Ï5|. The size of the condensate 
is limited by the initial number of partides in the gas, the 
initial size of the cut-off Eb(0), and the initial energy E 
(for 7 <C t~TÍ,). For 7 > t"~/k only few collisions will 
occur while the cut-off is ramped down. The number of 
partides that reach the condensate is therefore mainly 
determined by the collision rate. As can be seen from 



the Fig. there is a value for the ramp rate 7 which 
maximizes the number of partides transferred into the 
condensate per unit time, and therefore optimizes the 
cooling process under the assumption that additional loss 
rates from the trap do not change while the gas is cooled. 

We also performed some evaporative cooling simula- 
tions for a gas in a 3D-box. Because the density of states 
in that case does not rise as quickly as in the harmònic 
oscillator particle energies are changed more during a col- 
lision than in the harmònic oscillator. Therefore it is also 
possible to evaporatively cool a gas in a box quickly al- 
though the collision rate does not rise as much as in the 
harmònic trap. 

V. CONCLUSIONS 

We have simulated stationary and non-stationary prop- 
erties of a Bose gas in a trapping potential with finite 
number of partides in the framework of the Quantum 
Boltzmann Master Equation. 

For a gas confined in a 3D-box we have found that the 
number of partides in the condensate at a given temper- 
ature is larger than expected from the thermodynamic 
limit. We have also computed the mean collision time 
of partides in the gas. Comparison with the classical 
result shows that boson statistics tends to decreases the 
collision time close to the critical temperature. Calcu- 
lations of fluctuations in the number of partides in the 
one-particle ground state have shown that the Standard 
deviation increases almost linearly with temperature un- 
til the critical temperature is reached. For temperatures 
above T c the Standard deviation decreases again, and the 
distribution becomes Poissonian for high temperatures. 
We have also found that population is transferred at a 
time scale of the order of the collision time which is im- 
portant for the range of validity of the ergodic form of 
the QBME. 

Our simulations of a Bose gas in an isotropic harmònic 
trap were restricted to the ergodic form of the QBME. 
In contrast to the 3D box the number of partides in the 
condensate is decreased relative to the usual continuum 
(thermodynamic) limit at a given temperature. We found 
that the mean collision time decreases significantly as 
temperature reaches the critical point from above. This 
is due to the increase of the density, as soon as the ground 
state is macroscopically occupied. Simulations of evap- 
orative cooling have shown that there is an ramp rate 
to lower the cut-off energy of the trap with the goal of 
transferring as many partides as possible per unit time 
to the ground state. 

The present formalism is readily extended to includc 
mean field effects, and pumping and loss of partides from 
a degenerate Bose gas. This is relevant for modeling atom 
làsers based on collisions |Ï5],[Ï6) . 
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We insert expression (|A|) into the QBE @, divide by 
Ae, replace the notation (rq) by /(e^) and the ^ Ae 3 by 
/ de 2 des de^ and obtain finally 

, . , . 8nma 2 f , , , . . 
P(ei)/(ei) = — — Y- / de 2 de 3 de 4: ò(ei + e 2 - e 3 - e 4 ) 



ir 2 h 3 

{-/(ei)/(e 2 ) + /(e 3 )/(e 4 )} p(e min ). (A5) 



APPENDIX A: THE CLASSICAL LÍMIT 

To connect the present paper with Ref. we briefly 
rederive the classical Boltzmann equation with the er- 
godic approximation from the QBE (|Ï8|). We assume 
the distance between energy levels small compared to 
the mean energy of a particle so that the sum can be 
replaced by an integral. In the classical limit we get for 
the density of states at energy e 



P(e) 



J d 3 pd 3 x5^e~-U(x)-^j, 



(Al) 



where Í7(x) is the trapping potential. The degeneracy of 
the coarse grained one-particle states g(e) is connected 
to the density of states by g(e) — Aep(e). We replace 



£**(x)^(x') 



Ae 



(2tt^) 3 

i(x-x')pi/?i 



d 3 Pi 5[ei-U(x)-^ 



(A2) 



where e, = Tiují. The factor Ae in Eq. (A2) ensures the 
normalization of the sum over the wavefunctions to gifii). 



Inserting replacement Eq. ( |A2| ) into | f/1234 1 and integrat- 
ing over x' yields a ó-function of the four momenta times 
(27r7ï) 3 . Integrating over p 4 i.e. setting p 4 = P1+P2 — P3 
we obtain 



^ \Ui2m\ 2 5{e x + e 2 - e 3 - e 4 ) 



lel,2€2 
3e3,4e4 



16tt 2 7íV Ae 4 
m 2 {2-kTi)- 



d 3 pi d 3 p 2 d 3 p 3 



e l - C/(x) - |Q í(ei + e 2 - e 3 - e 4 ). (A3) 



We define the total momentum P = pi + P2 and the 
relative momenta q' = (pi — P2)/2 and q = (p 3 — p 4 )/2. 
Integrating over the azimuthal àngels of the two relative 
momenta q and q' and over the length of the relative 
momentum q and calculating the remaining integral sim- 
ilarly to M we obtain 



22 IC/1234I 2 5(ei + e 2 - e 3 - e 4 ) 



iei,2e2 
3e3,4eï 

2mAe 4 a 2 



n 2 h 2 



p(e m in)<Kei + e 2 - e 3 - e 4 ). 



(A4) 



In Eq. (A5) the (1 + /) factors are neglected by assuming 
that in the classical limit the mean occupation of a quan- 
tum level is much less than one. Setting a = 877a 2 this 
is the ergodic form of the classical Boltzmann equation 
from Ref. É. 
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FIG. 1. Condensate fraction versus temperature in thermal 
equilibrium for the 3D square well potential. (a) thermody- 
namic limit, (b) grand canonical solution for N = 500 (sòlid 
line) and results from the simulation (+), (c) grand canon- 
ical solution for N = 100 (sòlid line) and results from the 
simulation (o). 
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FIG. 2. Occupation numbers for the ground state N c and 
the first excited state n\ against the total number of partides 
TV in thermal equilibrium for T = 0.5 T c . 



FIG. 9. Distortion of an ergodic distribution into a non er- 
godic one for the 3D box. Particle distribution of the depleted 
levels (P(nf)) and of the filled levels (P(n{)) at time t after 
the distortion. Simulation for TV = 500 at T = 0.4T C . 



FIG. 3. Mean collision time per particle versus tempera- 
ture for TV = 500 of the 3D box without ergodic assumption. 
Result of the simulation (+) and result for í co n = (anvT)~ 
(dashed line). The time scale is normalized to {nav\2/ TV) -1 . 
i>i is the magnitude of velocity of a particle in a first excited 
state as defined in the text. 



FIG. 10. Condensate fraction versus temperature in ther- 
mal equilibrium for the 3D harmònic oscillator. (a) thermody- 
namic limit, (b) grand canonical solution for TV = 500 (sòlid 
line) and results from the simulation (+), (c) grand canon- 
ical solution for TV = 300 (sòlid line) and results from the 
simulation (o). 



FIG. 4. Mean collision time per particle versus tempera- 
ture for TV = 100 for the 3D box. Result of the simulation 
without (b) and with (a) the ergodic assumption, result for 
ícoli — (aftVT) -1 (dashed line). The time scale is normalized 
to {ncrv^l/ N)~ l . v\ is the magnitude of velocity of a particle 
in a first excited state as defined in the text. 



FIG. 11. Total energy of the system versus temperature. 
(a) data for harmònic oscillator, (b) data for the 3D box each 
with TV = 500. Energy is normalized to the level spacing Ae. 
(+) results from the microcanonical simulation; (sòlid line) 
result of the grand canonical calculation. 



FIG. 5. Probability distribution of partides in the conden- 
sate for the 3D harmònic box without the ergodic assumption. 
Results of the simulation (bars) and fits (sòlid lines) are cal- 
culated for TV = 500. 



FIG. 6. Probability distribution of partides in one of the 
first excited states for the 3D box without the ergodic assump- 
tion. Results of the simulation (bars) and fits (sòlid lines) are 
calculated for TV = 500. 



FIG. 7. Fluctuation of the condensate fraction versus tem- 
perature in thermal equilibrium for the 3D square well po- 
tential. Results of the simulation for TV = 500. The crosses 
give the results from the numerical summation of Eq. ([l2"|). 
The dashed line is y/Nl which would be equal to a(N c ) if the 
fluctuations in the condensate were Poissonian. 



FIG. 8. Buildup of the condensate for the 3D box for 
TV = 500. The energy is chosen such that in equilibrium 
T — 0.5T C . The time scale is normalized to (navi2/N)~ 
in thermal equilibrium. The dashed line is a fit of the form 
TV C (1 - exp(-í/r)), with r = 0.0013 and TV C = 368, as ex- 
plained in the text. 



FIG. 12. Mean collision time per particle versus temper- 
ature for TV = 500 for the harmònic oscillator. Result of 
the simulation (sòlid line) and result for í co n = (o-n^ur) -1 
(dashed line). The dotted line shows the collision time with 
the assumption of a fixed density equal to the ground state 
density. The time scale is normalized to (n < lavo2/N)~ 1 . vo 
is the amount of velocity of a particle in the ground state as 
defined in the text. 



FIG. 13. Total number of partides TV and number of par- 
tides in the condensate TV C for 7 = 1/10 (sòlid line) 7=1/2 
(dashed line) and 7 = 3/2 (dotted line) against time t. 7 
is the time constant from Eq. (BG) normalized to crn°«o2/TV. 
The time t is normalized to (a^Vo2/N)~ 1 . 



FIG. 14. Mean collision time í co u versus time t for 7 = 1/10 
(sòlid line) 7 = 1/2 (dashed line) and 7 = 3/2 (dotted 
line) against time t. 7 is the time constant from Eq. ( |3C| ) 
normalized to an1vo2/N. The time t is normalized to 
(an° h v 2/N)-\ 



FIG. 15. Size of condensate divided by time to reach 90% 
of the final size of the condensate versus time constant 7. 7 
is the time constant from Eq. (|3o|) normalized to an1vo2/N. 
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